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In a previous paper, we discussed estimation of the parameters of a 
single tone from a finite number of noisy discrete-time observations. In 
this paper, we extend the discussion to include several tones. The Cramer- 
Rao bounds are derived and their properties examined. Estimation 
algorithms are discussed and characterized. 

I. INTRODUCTION 

In a previous paper, 1 we reported on the estimation of the parameters 
of tones from a finite number of noisy, discrete-time observations and 
described the case of a single complex tone. In this report, we discuss 
the situation when the signal consists of several, say k, tones, either 
real or complex. By real signal we mean 

it 

S(t) = Y, °i COS W + 0<)« 

t"=l 

The corresponding complex signal is of the form 

s(t) + js(t) = E biexpljM + di)l, 

i—l 

where s(t) is the Hilbert transform of s(t). 

A computer observes, through the A-D converters, noisy versions 
of the signal, X(t), and possibly its Hilbert transform Y(t). That is, 
samples are taken of 

x(t) = s(0 + w(t), (i) 

and 

Y(t) = S(t) + W(t), (2) 

where W(t) and W(t) are the noise and its Hilbert transform, 
respectively. 
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The observations are made at times denoted t„. The computer will 
process one or both sample vectors : 

X = IX , X„ • • -, X N -{] T and Y = [To, Y u ~*t Y»-t] T , 

where T denotes matrix transpose, 

X n = X(t n ), (3) 

and 

Y n = Y(t R ). (4) 

We assume the noise samples, W n and W n , are independent, zero- 
mean, gaussian random variables with variance a 2 . 

Let a be the p-element vector of unknown signal parameters. We 
assume all signal parameters are unknown, so that p = 3k, and use the 
convention : 

OiZi-2 — Wit /r\ 

tt3t-i = hi, 
and 

a 3 i = Qi, i = 1 to k. 

This model describes several situations. The real signal may be 
received from a data set during a test or it could be a probe signal used 
to characterize a data-transmission channel. The real and imaginary- 
parts of the complex signal could occur as the result of in-phase and 
quadrature modulation processes, as described by Palmer. 2 The imagi- 
nary part of the complex signal could be the output of a 90-degree 
phase-shift network (Hilbert transformer) through which the real 
signal is passed before the sampling is done. This is done in certain 
types of data sets that use all-digital means to demodulate received 
signals. Samples of the complex signal are easier to process because of 
the absence of negative frequency components, as we show below. The 
model also applies to certain mathematically equivalent, phased-array 
radar problems, such as the one described in Refs. 3, 4, 5, and 6. 

There are two main aspects to the problem of estimating the pa- 
rameters of the signal: lower bounds to estimation accuracy and 
algorithms for doing the estimation. In the next section, the properties 
of the Cram&'-Rao (c-r) bounds are explored. There are many other 
bounds that could be applied but we have only examined the c-R 
bounds. Section III describes and evaluates some approximations to 
maximum-likelihood (ml) estimation. In Ref. 1, we found that when 
the signal consists of a single complex tone (fc = 1), then ml estimates 
can be obtained with any desired accuracy. When several tones are 
present, ml estimation is sufficiently complicated that suboptimum 
alternatives are attractive. 
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II. CRAMER-RAO BOUNDS 
2.1 General theory 

Maximum-likelihood estimates of signal parameters are unbiased at 
high signal-to-noise ratio (s/n). 4,7 We will develop estimation algo- 
rithms that have very little bias, so we have only studied the c-R 
bounds to unbiased estimation accuracy. Even when an estimator has 
some bias, the unbiased bounds serve as useful goals for estimation 
accuracy. Since the accuracy of ml estimates approaches the unbiased 
c-r bounds at high s/n, the unbiased bounds also show what could be 
done if exact ml estimation algorithms were used. 

We found in Ref. 1 that low s/n is that range of s/n where estima- 
tion anomalies occur. None of the known bounds seem to be very tight 
under these conditions. 

The first property of the c-r bounds that we consider is one for which 
we need the following general notation. 

Let V be a "signal" vector whose typical component is of the form 

V„ = £ &«<(«<, fc, n). (6) 

Notice that each </,•(•) has an associated level, 6„ and is a function only 
of n and the ith set of unknown parameters. Time does not necessarily 
enter into the Qi( m ) functions. Let X be a noisy observation of V. 
Assume the noise is additive, multivariate normal with zero mean and 
correlation matrix R _1 . If the noise vector is W, then 

X = V + W, (7) 

and the probability density function of X given V is 

/(X/V) = ^L exp[-HX - V)'R(X - V)], (8) 

where N is the dimension of V and the T denotes transpose (see Ref. 8, 
page 207). 

The c-r bounds require certain regularity conditions on V, which 
are satisfied by our model. 9 The bounds are the diagonal elements of 
the inverse of the Fisher information matrix, J, whose typical element 10 
is: 



Jab — —E 



<V log/1, ( 9 ) 



dabda a 
where E { • } denotes expected value of { • } . The bounds are : 

Var {d a - a a ] ^ J«*, (10) 
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where «/°° is the ath diagonal element of J -1 and & a is an unbiased 
estimate of a a . 
It is easy to show 11 that 

dab oota 

We now present a few theorems that characterize the c-r bounds. 
We assume J is not singular, for reasons discussed below. 

Theorem 1 : The C-R bounds to unbiased estimation of the parameters o>,- 
and 0, of V are functions of b, but are independent of the other levels, 
bj,- j y£ i. The bound to unbiased estimation of a level, bi, is independent 
of all the levels. 

Proof: Equation (11) is equivalent to 

^ 6 =£Lfln m ^^, (12) 

n m da b aa a 

where R nm is an element of R. 

The elements of J that are functions of the parameters of g, and gj, 
using the convention given by (5) and the notation g,(n) = gi(ui, 0., n), 
are: 

Js.-2,3i-2 = bibj 2- E «n„ -3 « (laa; 

Jii-w-i = 6.11 R nm Bj £^ ft(m). (13b) 



71 m 



t n,^vi? fain) dgjjm) m , 

Jsi-i.tj = bib, L E Rnm du dd • Uac; 



n m 



/ikw = h E E B.^(n) %^ • (13d) 

/m-i.V-i = E E Rnmgi(n)g } (m). (13e) 

n m 

/,<-!.., =b,II B.^*(n) ^^ • (13f) 

Jii.ti-i = bibj 2, E «»-» — ^ a -- ' "*® 



n m 



/«,W =biZZ Rnm ^SM gj (m). (13h) 



n m 



j£«(w) gg/(m) n o:x 

•/M.V = OiOy 2, 2- Rnm —^ QQ~ ' ^ ' 



n m 
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An examination of (13a) through (13i) shows that the submatrix of 
J has the form D.Q.yDy, where 



D,= 



\bi 





0] 





1 











bi 



(14) 



and the matrix Qy is not a function of any &,. It follows that J has the 
form 



where 



and 



D = 



Q = 



J = DQD, 




p X 


1 


D 2 




" 


D*. 



Qi* 



is a matrix whose elements are all zeros. From (15), 

J-i = D-iQ-iD" 1 , 

from which the theorem follows. For example, 
Var {a x - an} ^ Q n /bl 



(15) 



(16) 



(17) 



(18) 



(19) 



This theorem is not entirely new. It is alluded to in Ref. 6. However, 
this form of the theorem shows that, contrary to Ref. 6 and popular 
opinion, precisely known sampling times (or antenna element spacing 
in the equivalent radar problem) are not necessary for the theorem to 
hold. 

The theorem is true whether or not the noise samples are independent 
and regardless of the sampling times. Of course, if the sampling times 
are not known, then the c-r bounds cannot be accurately calculated, 
but that does not obviate the theorem. 

It should also be clear that the number of unknown parameters is 
unimportant to the theorem. Clearly the theorem holds if, for example, 



and if 



<7,-(wi, 6i, n) = cos (u)it„ + Bi), 



t n = nT. 



Theorem 2 : The bounds associated luith the parameters of the first k tones, 
when there are k + m tones, are not less than the bounds when there are 
only k tones. 
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Proof: The matrix J is always positive semidefinite. Thus, if it is not 
singular, it is positive definite. 

Suppose J is the Fisher information matrix for k + m tones and is 
partitioned so that J* is the J matrix associated with k of the tones. 
This partitioning is always possible. Then write 



J = 



h | K ' 



(20) 



Since J is positive definite, so are J k and J m . 
Write the inverse of J in the form 

W 



W r 



(21) 



where J* and U are both 3A; by 3A; matrices. Theorem 2 is true if 

U ^ Jf 1 , (22) 

which means U — Jf * is positive semidefinite and which we now prove. 
Using the fact that 

JJ- 1 = I, (23) 

one can show that 

U = [J* - KJ" 1 ^]- 1 . (24) 

Observe that KJ^K 7 " is positive semidefinite. That is, 

SJ;V ^ 0. (25) 

Since J fc and hence U are positive definite, (24) and (25) imply that 
U" 1 S Jk, which implies (22). 

Another implication of the proof of Theorem 2 is that the bounds 
for p of p + m unknown parameters are not less than the bounds 
when only the first p parameters are unknown. 

This theorem is also not entirely new, although we have not seen it 
stated before. A restricted version of the theorem is mentioned in 
Ref. 12, (page 33), and Problem 2.4.23 in Ref. 10 hints at this kind of 
result. 

The theorem depends upon J being nonsingular. It is easy, but 
tedious, to show that if the signal vector is composed of samples of the 
real or complex signal described in the introduction and only two tones 
are involved, then J is singular only if the two tone frequencies are 
equal, modulo 2ir/T, where T is the intersample time. (Remember that 
a real tone has a component at +0, and another at -co,-.) 
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We have not been able to prove this result for an arbitrary number 
of tones, but all of our calculations of various J matrices support the 
hypothesis that J is singular only if two or more of the tone frequencies 
are equal, modulo 2t/T (assuming N is large enough). 

When two of the tone frequencies are equal, the receiver is receiving 
one less tone than expected. In this paper, we assume that the correct 
number of tones, k, is known and that all of the frequencies are distinct. 

2.2 Equally spaced samples and independent noise 

We now concentrate on the problem described in the introduction. 
Assume all noise samples are independent with variance a 1 . That is, 



R = -,I, 



where I is an identity matrix. 
Define 



(26) 

(27) 
(28) 

(29) 

As is mentioned in Ref. 1, the time of the first sample, t , has an effect 
upon bounds and estimation accuracy. We have ignored that problem 
in this paper and taken to to be zero. 
The signal vector is 



and 



Hn = Z bi COS (Uitn + 6i), 



Vn = £ bi Sin (tOitn + Bi), 

i-1 



tn = nT; n - 0, 1, •■■,#- 1. 



V n = 



V n -N 



n = to A r - 1 
n = N to 2N - 1 



or 



7» = /»»; n = 0toN-l 
Then a typical element of J is 

- L N T [fyndMn dVnd Vn 1 

" a 2 „ = o |_ da a dab da a dab J 

The v n terms are dropped if the signal is real. 
Let M («, 6) be the matrix defined by 



(complex signal) (30) 
(real signal). (31) 

(complex signal). (32) 



M(w,0) = 



T 2 £ n % cos An - T £ n sin A„ T E n cos A„ 
T^n sin A n L cos A„ £ sin A„ 

7 1 E n cos A„ — E sin A„ E cos A„ 



(33) 
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where 

A n = nuT + d; n - to N - 1. 
Let P,-y be the matrix defined by 

P„ = M(oJi - a}j,di - $f) 
and let P be the p-by-p matrix defined by 

Pii • ■ • Pit 



P = 



'hi 



Let 



B = 



1 
0-10 
1 



(34) 
(35) 

(36) 
(37) 



Define a matrix Q„ by 

Qa = i[M( Wi - ^, di - 6j) - M(o>,- + wy, Bi + y )B] (38) 
and a matrix Q by 

"Qn • • • Qi» 



Q = 

Q*i • • Qkk 

Then it can be shown (13) that J is given by : 



or 



J-^DPD 



J = pDQD 



complex tones 



real tones. 



(39) 



(40) 



(41) 



Theorem 3: When the signal consists of two equal-level complex tones, 
the C-R bounds for the same parameters (e.g., the two frequencies) are 
equal. In other words, the mutual interference is reciprocal. 

Proof: The J matrix is 

because P„ = P 22 = M(0, 0). Observe that Pf 2 = BP 12 B because 
M r (oj, 6) = BM(co, 0)B. Thus, J -1 has the form 



J ~ ff [ D^JLw^ V JL Dr\T 



(43) 
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where 
Hence, 



V = BUB. (44) 



Uij\ = \Vtfl (45) 



and Uu = Vu, which proves the theorem. 

Theorem 4: The bounds for two tones, real or complex, are periodic in 
0! and 6 2 with period ir. 

Proof: The theorem follows from the easily checked fact that 
M(w, 9 + v) = -M(a>, 0). 

Theorem 5: The bounds for real or complex tones are periodic in each 
frequency with period 2w/T. 

Proof: The theorem follows from the fact that M(o» + 2n-/T, 9) 
= M(w, 6). 

Theorem 6: The bounds associated with complex tones depend upon the 
difference frequencies and phases but not upon the absolute values. 

Proof: The theorem follows from (35), (36), and (40). 

It is, in general, tedious to invert J and obtain formulas for the 
bounds. However, it is a simple matter to have a computer calculate 
the elements of J and its inverse. We have done this to obtain a better 
understanding of the bounds. 

A number of illustrative curves are given in Ref. 13. In the interest 
of brevity, we will present only two of the figures here. 

The main thing we learned from the calculations is that there is a 
critical frequency separation, 4lt/NT, associated with multitone c-R 
bounds. In Ref. 1, it is shown that when a single complex tone is 
present, the bounds are independent of the frequency of the tone. When 
more than one complex tone is present, the bounds approach the single- 
complex-tone bounds when the minimum frequency separation (modulo 
2v/T) exceeds the critical frequency. The multitone bounds increase 
rapidly as the minimum frequency separation goes below this critical 
frequency. 

This rule applies to a single real tone if it is considered to be two 
complex tones, one at a frequency, say, of wi, and one at — u\. Thus, 
if the frequency of a single real tone is less than 2t/ NT, modulo ir/T, 
then its c-R bounds are much larger than the corresponding single- 
complex-tone bounds. 

In all cases, the multitone bounds depend upon the tone phases, as 
might be expected. 
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Fig. 1 — Frequency estimation bound vs A/ for center of three equally spaced 
complex tones with worst relative phase and 20 dB s/n. N is 128. \/T is 4000 Hz. 
Corresponding single and double tone bounds also shown. 

Figure 1 illustrates the critical frequency for frequency-estimation 
bounds. The worst phase, i.e., the phase that gives the largest frequency 
estimation bound, was used at each difference frequency. Figure 2 
shows the critical frequency effect upon the frequency-estimation 
bound for a single real tone. 

To facilitate comparisons, in all figures we used a sampling frequency 
of 4000 Hz for complex tones and 8000 Hz for real tones. Thus, in both 
cases, the unknown tone frequencies are assumed to fall in the range 
of to 4000 Hz. 

III. ESTIMATION ALGORITHMS 

3.1 General 

The ml estimation procedure is conceptually simple. Given that a 
sample vector, X, is received, the ml estimate of the parameter vector, 
a, is the value of a that maximizes the p.d.f. of X. That is, a maximizes 
/(X/V). a may not be unique. Maximum likelihood estimation of the 
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parameters of a single complex tone has been shown to be relatively 
easy to implement. 1 It was shown in Ref. 1 that single complex-tone 
ml estimators have variances almost equal to the or bounds over a 
wide range of s/n. No other unbiased estimators could do significantly 
better over that range of s/n. 

Maximum likelihood estimation when several tones are present is 
much more difficult to implement. However, we show below ways to 
approximate ml estimation. We start the discussion with complex 
tones and examine a practical approximation to ml estimation, the 
resulting bias effects, the use of window functions to reduce bias, and 
a time-saving interpolation algorithm. Then we briefly discuss how the 
ideas and results apply to real tones. 

Recall from Ref. 1 that we seek to maximize the function 

L = | E (XnM- + Y nVn ) - i Z ( M 2 „ + vl), (46) 

where X n and Y„ are as defined in (3) and (4). 

After carefully arranging the terms, we obtain the likelihood func- 




300 400 
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700 



Fig. 2 — Frequency estimation bounds vs frequency for single real tone at 20 dB 
s/n and worst phase. 
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tion in a form similar to that used in Ref . 1 : 
L= E { 26, Re |V*'A (»<)]-&?} 



- h L £ &<&m £ cos (na>,T - na n T + 0,- - m ), (47) 



where 



A(«) = i £ (X. +jT„)e-'- r . (48) 

-<V „=o 

L as given by (47) has two main terms and would be difficult to 
maximize by a simple program. It can be done, but a lot of work is 
involved. We notice, however, that when there is only one tone (k = 1), 
the second term of (47) vanishes. Also, when N is large and k > 1, the 
magnitude of the second term is still relatively small and does not in- 
volve the data. Thus, we are led to drop the second term in L and 
maximize the remainder. This, of course, will only give ml estimates 
when A; = 1 and will give "almost ml" estimates otherwise. 

3.2 An almost ML algorithm 

Suppose the cross-product terms in (47) are dropped. Then to make 
estimates, we need to maximize 

U - £ 2bi Re \jr**AM] - b\. (49) 

From Ref. 1, each frequency estimate, <$, maximizes | A (u>) \ . Then the 
corresponding level and phase estimates are 

S<= I A (4) I (50) 

and 

h = arg [A (»*)]. (51) 

The function | A (w) | has many maxima and large peaks near the 
frequency of each tone. Thus, the frequencies of these large peaks, as 
illustrated in Fig. 3, are taken to be the frequency estimates, &,. Due 
to the periodicity of | A (a>) | , all the a>, should be confined to a range 
no wider than u, = 2-tr/T to avoid ambiguous frequency estimates. 
Normally the range (0, 2r/T) is used. When real tones are involved, 
the range should not exceed t/T. 

3.3 Bias 

Consider the case of only two tones. An example of | A (w) | when the 
noise power is zero is shown on Fig. 3. 
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FREQUENCY (oj) 

Fig. 3 — Shape of \A («) | from two complex tones of equal phase and level, without 
noise. N is 16. 

The figure has large peaks near a>i and o>2. The peaks in the example 
are actually both displaced away from the average of the two fre- 
quencies. Thus, the penalty for neglecting the cross-product term in 
(47) is a bias in estimates of frequencies and levels. 

The frequency and level bias in the zero-noise case is easily cal- 
culated. An example of such calculations is shown on Fig. 4. The figure 
shows the dependence of frequency estimation bias on the difference 
frequency (A/) of the two tones. When two tones have almost the same 
frequency, the two large peaks merge into one at a frequency equal to 
the average of the two tone frequencies. This accounts for the negative 
slope of — \ at low A/ on Fig. 4. There is also a dependence upon the 
difference phase (A0). 

Figure 4 shows the bias for one of the two complex tones. The bias 
for the other has the same magnitude but opposite sign. In general, 
the magnitudes of the biases for two tones are not equal. However, 
they are equal when the two tones are equal-level complex tones. 

3.4 Window functions 

In discrete Fourier transform (dft) work, window functions (also 
called weighting functions) are often used to minimize the effects of 
one tone upon another. The modification of the dft of samples of one 
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Fig. 4 — Bias in peak frequency of |A(w)| from two equal-level, equal-phase com- 
plex tones vs difference frequency, without noise. N is 64. 

tone by the presence of samples of another tone is called leakage. See 
Rife and Vincent 14 for a discussion of leakage and how window function 
will reduce it. 

In the time domain, a window function (or time window), say h(t), 
is characterized by its samples, h(nT). In use, each data sample, 
X n +jY n = Z„, is multiplied by h„ = h(nT) before A (w) is computed. 
Thus, A (w) becomes 



A(C0) = -^ L /inZne-""" 3 '. 



(52) 



When a window function is used, the bias in the frequencies of the peaks 
of | A (w) | is modified. If a good window is used, the bias can be greatly 
reduced. The penalty, as we see below, is an increase in the variance 
of &i and 6,. Palmer also reported this penalty in Ref. 2. 

In the context of the dft, window functions can be written in the 
form 



M 



h n = 1 + £ di cos (2Tin/N). 
i-i 



(53) 



The number M, which can be assumed to be less than N /2, and the di 
define particular windows. With h n in this form, 

I N-l 

Jj L hn = 1. 

JV n = 
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Table I — Values of d { in ascending order of / for 
various windows 



Hanning 


Standard 


Taylor 


-1 


-1.43596 


-1.03538 




0.497536 


0.0824936 




-0.061576 


-0.00116197 

-0.00188862 

-0.00123387 

-0.000671595 

-0.000275885 



A window that is better than many at reducing bias is the one 
identified by Rife and Vincent 14 as g»(t). We call this the standard 
window. Another useful window is one of the Taylor windows. 14 These 
windows are defined in Table I. 

Figure 5 is an attempt to summarize the way window functions 
affect bias. The curves on the figure compare upper bounds to the bias 
associated with each of the previously defined window functions. The 
curves were obtained by computing at each frequency the bias at the 
worst phase (the phase that gave the largest bias) . The resulting curves 
were flattened as indicated for the Taylor curve. 

Figure 5 shows the Taylor window does the best job when the tone 
frequency separation is small. At large separations, however, the 
standard window does much better. The figure also shows how bad 
the bias is if no window is used. 

Windowing increases the variance of frequency and level estimates. 
It can be shown 13 that the increase in variance is related to the function. 



1 JV-1 

— — y f,2 

JN n=0 



(54) 



It is easy to show that 



= 1 + i £ dl 



if 



N > 2M. 



(55) 



Thus, 7] is not a function of N. Simulations verify that larger rms 
errors are associated with larger values of tj. Some values of rj are 
tabulated below. 

Window 7j 



None 


1.00 


Hanning 


1.50 


Taylor 


1.54 


Standard 


2.16 
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Durrani et al. call the sum r) a dispersion factor. 15 They have com- 
pared many windows and have tabulated their parameters, including 
dispersion factors. Other windows are mentioned in Blackman and 
Tukey. 16 

The data on Fig. 6 shows the general effects of windows on rms 
errors when a single complex tone is present. The Hanning window 
produces almost the same rms error as the Taylor window and is not 
shown on the figure. 

Bias contributes to rms errors more than variance does at high s/n, 
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Fig. 5 — Magnitude of frequency estimation bias for two equal-level complex tones 
using window functions. Curves are leveled as described in the text. Worst-phase was 
used at each frequency. N is 64. 
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Fig. 6 — Simulation results showing the effect of window functions on estimation 
variance with a single complex tone. N is 64. 

while estimation variance controls rms errors at low s/n. Thus, while 
a given window may produce lower rms errors than another at high 
s/n, the roles may be reversed at low s/n. The "best" window for a 
given application will, therefore, depend upon the tone frequency 
spacings, the expected s/n, and possibly other factors. Figure 7 illus- 
trates this point. On the figure, the Taylor window is best at 10 dB 
s/n, but the standard window is best at 40 dB s/n, where the bias 
associated with the Taylor window causes the rms error curve to level 
off. 

3.5 Interpolation 

Maximization of \A («) | involves a search routine. A two-step 
algorithm that has a coarse search and a fine search was described in 
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Fig. 7— Simulation results showing combined effects of bias and variance on 
frequency estimation for one of two complex tones. Frequency difference is 380 Hz ; 
worst-phase was used for each window. N is 64. 

Ref. 1. Fine searches are time consuming. This can be serious if com- 
puter time is important. One way to trade accuracy for speed is to use 
an interpolation algorithm on the dft of the input data to arrive at 
frequency and level estimates. 

Rife and Vincent developed several interpolation algorithms. 14 The 
one we investigate here is the following. 

Assume the output of the fft is the set : 



JS n = 



k m tO N - 1. 



(56) 



Suppose a coarse search is conducted over < k < N. This results in 
locating \A t \ which is the largest |A*| in the interval. Choose a = ±1 
such that | A i +a \ ^ | Au-a I 
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Let 
and 



ai = \Ai\ 



a 2 = \A 



l+a 



Assume the sampling frequency is 03, = 2r/T. 
The formulas from Rife and Vincent are : 



and 



6 = 



6 - I (I + aS) 



2tt(hX 



sin (ttX) \ 



M "1 ' 

1 + £ d n 8>/(8* - n 2 ) 

n = l J 



where the d n define a window and 

da 2 — Cid\ 



8 = 



C3CI2 + a x 



(57) 
(58) 

(59) 
(60) 

(61) 



The numbers C\, C2, and C 3 are given by Rife and Vincent in Table 
II for several windows. 

The interpolation formulas give estimates that are only a little 
worse than the fine search gives, rms frequency errors are typically 
increased by about 30 percent when interpolation is used, rms level 
errors increase less. 

When many tones are present, window functions can provide a 
satisfactory reduction of leakage as long as the minimum frequency 
separation is no less than about 8ir/NT. The data on Fig. 8 illustrate 
this point. The tone phases were all made random for these simulations. 
Thus, the points indicate the rms errors one might encounter in a 
working system. The bound shown on the figure is the (unbiased) c-r 
bound maximized over the possible phases of the center tone. 

We consider a real-tone estimation system to be equivalent to a 
complex-tone system if the two systems have the same useful band- 
width and the same frequency resolution. This means (i) the real 
sampling frequency is twice the complex sampling frequency and (it) 
the total sampling time, nt, is the same for the real tones as for the 



Table II- 


-Constants for computing delta 


in eq. (61) 


Window 


c, 


c 3 


c, 


None 
Hanning 
Taylor 
Standard 


1 
2 

1.96339 
3.6020 




1 

1.01643 

2.5862 


1 
1 

0.893534 
1.0317 
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Fi gt g_ Simulation results showing effects of Taylor window and interpolation 
algorithm upon frequency estimates of center tone of three equally spaced real tones. 
Center tone frequency is 2000 Hz. All have random phase. N is 128 and s/n is 20 dB. 

complex. For example, a real-tone system using 1/T = 8000 Hz and 
N = 32 is equivalent to a complex-tone system using 1/T = 4000 Hz 

and N = 16. 

The estimation algorithms described above for complex tones can 
be applied to real tones whose frequencies, in Hz, are in the range 
(1/iVT, 1/2T — 1/NT). The resulting accuracies are about the same 
as in the equivalent complex case. 

IV. CONCLUSIONS 

We have studied the problem of estimating the parameters, such as 
level and frequency, of several sinusoidal signals from a number of 
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noisy observations, taken at discrete-time instants. Gaussian noise and 
ideal analog-to-digital conversion were assumed. The nature of the 
problem led us to study the generalized Cramer-Rao lower bounds to 
estimation accuracy and maximum likelihood estimation. The com- 
plexity of maximum likelihood estimation algorithms led us to examine 
several algorithms that yield estimates that are almost, but not 
exactly, maximum likelihood estimates of the signal parameters. 

We were able to obtain estimators that have negligible bias, at least 
at high s/n. Thus, we considered in detail only the generalized c-r 
bound for unbiased estimators. Even when the resulting numbers are 
not, strictly speaking, lower bounds (e.g., when an estimator is biased), 
the unbiased estimation bounds can be considered to be desirable 
objectives for estimators. 

Several properties of the bounds were derived from the properties 
of the J matrix. Other properties, such as the existence of critical 
frequencies, were revealed from computations. 

The J matrix in the real tone cases is more complicated than in the 
complex cases and does not have quite the same structural properties. 
Thus, for example, the lower bounds for a single complex tone are not 
also lower bounds for the equivalent single real tone. On the other 
hand, the bounds for the case of many real tones approach the bounds 
for the equivalent complex cases when none of the real tones have 
frequency differences less than 2/ NT, modulo 1/27 7 (in Hz). 

The cases of many complex tones and of real tones present some 
difficulties. Maximum likelihood estimation is difficult to implement 
because of the presence of cross-product terms. To properly implement 
ml estimation, multidimensional search procedures over a nonconvex 
function would be necessary. We found that when the tone frequencies 
are separated far enough, the cross-product terms could be neglected, 
thereby permitting the use of a simple algorithm whose estimates are 
almost equal to ml estimates. 

The penalty for dropping the cross-product terms is a bias in fre- 
quency and level estimates. We found that the use of a suitable window 
function will reduce the bias to the point where it can be neglected 
when the minimum frequency separation of the tones is 4/ NT. Three 
window functions were discussed and compared. 

We found that the use of a window to reduce bias increased the 
variance of the frequency and level estimates. The rms error of fre- 
quency estimates is increased by about 35 percent with Taylor window 
and by over 100 percent with standard window. The use of the in- 
terpolation formulas increases rms frequency errors by another 30 per- 
cent or so. Level estimates are affected less by windows and interpola- 
tion. All of these figures apply when the s/n is above threshold. 
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